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Abstract 

In this paper, we present a method for the computation of the low-lying real eigenvalues of 
the Wilson-Dirac operator based on the Arnoldi algorithm. These eigenvalues contain information 
about several observables. We used them to calculate the sign of the fermion determinant in one- 
flavor QCD and the sign of the Pfaffian in = 1 super Yang-Mills theory. The method is based 
on polynomial transformations of the Wilson-Dirac operator, leading to considerable improvements 
of the computation of eigenvalues. We introduce an iterative procedure for the construction of the 
polynomials and demonstrate the improvement in the efficiency of the computation. In general, the 
method can be applied to operators with a symmetric and bounded eigenspectrum. 
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I. INTRODUCTION 



The Wilson-Dirac operator D\y, which is used in many recent lattice simulations to 
represent the fermionic part of the discretized action, has the following form 

4 

n—^,m ■ (1) 

Here n, m denote points in a four-dimensional hypercubical space-time lattice, a, (3 are Dirac 
indices, /i = 1,2,3,4 labels the positive directions and 7^ are the Dirac matrices. The hop- 
ping parameter k is related to the bare fermion mass; in particular /t increases for decreasing 
fermion masses. The link variables are associated with the links connecting neighbor- 

ing lattice points and represent the gauge field. In our investigations [H [2] the gauge field 
was in the fundamental representation of SU(3) for QCD with one quark flavor (one-flavor 
QCD) and in the adjoint representation of SU(2) (real 3x3 matrices) for supersymmetric 
Yang-Mills theory. The method presented here is, however, not restricted to a specific gauge 
group and can be applied also to other fermion operators. 

In the free theory, the eigenspectrum of Dy/ can be decomposed into a physical branch, 
consisting of the smallest eigenvalues, and the doublers, which become irrelevant in the con- 
tinuum limit [3j. Such a clear distinction of relevant and irrelevant parts is not possible in 
the interacting case. However, the lowest part of the spectrum still contains the most im- 
portant information. The low eigenvalue part plays a crucial role in spectral decompositions 
of the fermionic observables and the lowest eigenmodes allow for an acceleration of the 
inversion by deflation [5]. 

For several investigations, the Hermitian operator 75-Dvi^ can be used instead of Dw- The 
corresponding eigenvalue problem can also be solved with other iterative methods, but for 
the non-normal operator Dw the (restarted) Arnoldi algorithm [6] seems to be the optimal 
choiceE] 

The importance of the lowest eigenmodes of Dw has been the subject of several recent 
investigations, e. g. in |Hl [9]. Furthermore, their implications on the topology of gauge fields 
has been studied (e. g. in [TU| [TT]). even though Dw does not allow for a realization of chiral 
symmetry on the lattice. 

In several cases, numerical simulations of field theories with dynamical fermions require a 
reweighting of the observables with the sign of the determinant of Dw or of its Pfaffian. This 
sign can be obtained from the number of negative real eigenmodes jHllS]. The computation 
of the reweighting for one-flavor QCD and M = 1 super Yang-Mills theory was the main 
purpose of our investigations of the spectrum of Dw- 

On small lattices, the complete set of eigenvalues is accessible (see e.g. Fig. [T]). In a 
more realistic setup, strategies focusing the computation on the relevant small eigenvalues 
and accelerating the convergence are required. For lattice QCD, a polynomial approach 
focusing on the low eigenmodes of Dw has been presented in [13]. Within a mathematical 
framework, other methods based on polynomial transformations have been developed for the 
computation of a particular sector of a general eigenspectrum [T^ [T5] . We explain here a 
new strategy to obtain the lowest real eigenmodes of the Wilson-Dirac operator and show its 
impact on the efficiency of the computation. Our strategy allowed us to obtain the relevant 
part of the spectrum on lattices up to a size of 32^ x 64 lattice points. 

^ See, e.g., [7 for a detailed discussion of the effects of the non-normahty. 
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This paper is organized as follows. In the next section we explain the basic idea of the 



polynomial transformations. In Section III we present our specific strategy to obtain the 



necessary polynomial. Section IV contains some results and Section |V] a comparison with 
other methods. Further mathematical explanations and some practical considerations can 
be found in Section IVII 



II. ACCELERATION AND FOCUSING OF THE ARNOLDI ALGORITHM 
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FIG. 1. As an illustration of the method and the form of the spectrum to be expected, this 
figure shows all eigenvalues of the Wilson-Dirac operator Dw from 29 independent thermalized 
configurations of a dynamical simulation of the A/" = 1 super- Yang-Mills theory [2] on a 6^ X 8 
lattice (/3 = 1.6; k = 0.1575). The red points correspond to the eigenvalues computed with the 
peeling method described below. The gray background shows the region of the eigenvalues in the 
free theory with the same value of k. 



Let ^i)^ be the region that contains all eigenvalues of the operator. We are only interested 

( f) ( f) 

in a subset of eigenvalues enclosed in a region denoted by ^j^w' case, ^j^w chosen 

to be a prolate region surrounding all real eigenvalues smaller (or larger) than a certain value 
(e. g. fi^^ = {Ail |5> [Xi] \ <e,^ [Xi] > x^m} with e small). 

Fig. [T] shows that the spectrum of Dw contains large regions with a high eigenvalue 

density and a nonzero imaginary part. For an efficient calculation of the real eigenvalues, 

( f) 

it is crucial to exclude these regions and focus the computation on the eigenvalues in ^Y)^- 
The Arnoldi algorithm computes the eigenvalues starting from those with largest real part. 
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It calculates the eigenvalues in the region f^^^ = {A|3ft [A] > a;!?^''}, where xif^"* depends 
on the parameters of the algorithm and the eigenvalue distribution of Dw- 

Hence, a direct computation does not focus efficiently enough on f^/j^ since a lot of 
unwanted eigenvalues are calculated. However, an appropriate polynomial transformation 
Dw ~^ p{Dw) leads to a better overlap of with fi^"^]-,^^. The computation gets 

focused on the relevant part and a smaller number of unwanted eigenvalues is computed. 
The eigenvalues of Dw can be obtained from the eigenvalues or eigenvectors of p{Dw)- 

The second advantage of the polynomial transformation is an acceleration of the Arnoldi 
computation. The computation of an eigenvalue A converges faster, if this eigenvalue is 
better separated from the rest of the spectrum (compared to some average distance of the 
eigenvalues). Therefore, a polynomial minimized on Q^'^'^ \ {A} (with p{\i) fixed) leads to an 
acceleration of the computation of Aj (for details cf. [TT, 16]). An analytic solution for the 
absolute minimum on a general ^l^^'^ \ {A} is not available, but Chebyshev pjj and Faber 
polynomials |il5j provide approximate solutions of it. 

Since the algorithm starts from a random initial vector, it can happen that some eigen- 
values within f^^^ are not found in the Arnoldi iteration. Especially, some in a set of closely 
lying or exactly degenerate eigenvalues might be missing. This effect is considerably reduced 
by the polynomial transformation. 

For an appropriate polynomial, the focusing effect and the acceleration by far compensate 
the costs of the additional multiplications. Eigenvalues in the original spectrum obtained 
with a polynomial transformation are shown in Fig. [Tj 



III. THE PEELING TRANSFORMATION 

In previous investigations of the lowest real eigenvalues of D-w in lattice QCD 
certain set of simple polynomials has been proposed. It consists of power transformations 
of the form 

PpoweriDw; n,a) = {Dw + crl)"" , withn G N, a G R. (2) 
It has been shown to considerably improve the performance of Arnoldi algorithm. The 



effect of this transformation on a test eigenspectrum is illustrated in Fig. |2(b) The region 

of computed eigenvalues in the original spectrum gets a wedge like shape. Hence, the 

( f) 

computation is better focused on ^jj^- However, at larger n the focusing effect saturates. 

Based on these observations, we propose here the "peeling transformation" as an iteration 
of the power transformation. It consists in the following steps: 

1. The starting point is a power transformation with an additional renormalization factor 
ro G R, po{Dw]no,ao,ro) = ppower{Dw/ro]no,ao). 



2. For the resulting eigenspectrum, a new power transformation is chosen for a further 
focusing on f^fj^^,, 

Pi{Dw] no, (To, ro, rii, ai, ri) = ppower {po{Dw; tlq, do, ro)/ri; rii, ai)). 

3. This procedure is iterated until the polynomial Pat is obtained. 



A choice of parameters n and a is explained in Sec. VI 
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The effect of the further iterations on a test eigenspectrum is shown in Fig. [2j Clearly, 
the eigenvalues in ^^^l^ are made accessible by the transformation, while ^^^^ \ is 
compressed in a region close to the transformed zero. The polynomial resulting from the 
iteration is 

PN{Dw;no,ao,ro, . . .,nN,o-N,rN) = (■ ■ ■ {{Dw/ro + dol)"" + . . .)/rN + ctnT'^ , (3) 
with the free parameter n^, (Jj, and rjj^ The question of an optimal choice for these param- 



eters depends on the form of the eigenspectrum and is addressed in Sec. VI 



IV. REAL EIGENVALUES AND DETERMINANT SIGNS 

One of the goals of our calculations were the determinant signs for numerical simulations 
of one-flavor QCD and Pfaffian signs for the supersymmetric Yang-Mills theory. In order to 
realize small pion or gluino masses, both theories were simulated within a parameter regime 
where very small and negative eigenvalues appear. Except for the real modes, all eigenvalues 
of D\y appear in complex conjugate pairs. Thus the determinant and Pfaffian signs depend 
only on the real negative eigenvalues; in particular 

det(Diy) = n I ^ I' n (4) 

{AeC|Q[A]>0} S[A]=0 

see [TTl ITS] for more details. For one-flavor QCD, Fig. |3] illustrates the distribution of the 
lowest eigenvalues for two different k values. 

This method, based on a direct computation of the real negative eigenvalues of Dw, 
turns out to be more efficient than the previously considered "eigenfiow" approach, where 
determination of the sign is based on the eigenvalues of ^^D^. This Hermitian matrix allows 
to compute its eigenvalues by means of simpler computational methods. However, in order 
to obtain the determinant (Pfaffian) sign the eigenvalues have to be computed at several 
different k values HZlllSjg 

Depending on the parameters of the simulation, in particular on the value of the hopping 
parameter k, we obtained up to 50% negative signs in the simulation of one-flavor QCD. 
With increasing k, approaching its critical value Hc, the number of negative signs increases. 
For values of k, which have been used for measurements of the particle spectrum, however, 
the number of negative signs was well below 10%. 



V. COMPARISON OF THE POLYNOMIAL TRANSFORMATIONS 

To demonstrate the importance of the different steps in the peeling method, we compare 
the performance of several different peeling polynomials with power polynomials. The com- 
putation time needed to obtain a number of wanted eigenvalues allows for a simple and clear 
representation of this performance. In the present case, n'^' is the region of all eigenvalues 
with an imaginary part whose absolute value is smaller than 0.05. The polynomials were 
constructed as described in Section IVIl 



^ The parameter can be absorbed into a redefinition of the (Ji and an overall rescaling. 
^ This method is similar to the computation of all real eigenvalues in |19| . 
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FIG. 2. The figures on the left hand side show the result of the polynomial transformation applied 
to a test spectrum of equidistant eigenvalues that fill a rectangular region. The small region colored 
in light blue corresponds to ^£)^ ■ The figures on the right display the parts in the original spectrum 
computed successively in the polynomial Arnoldi iterations. In the first figures no polynomial is 
applied. The polynomials applied in the other figures are listed below each of them. 
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FIG. 3. Distributions of the lowest real eigenvalues of the even-odd preconditioned Wilson-Dirac 
operator in simulations of one- flavor QCD. When k is increased, the distribution is shifted and a 
larger peak in the negative sector is observed. 



The eigenvalue computations depend on the number of computed eigenvalues, the size 
of the available eigenspace, and the maximal number of iterations. We have varied all of 
these parameters for the comparison of the polynomials. Fig. |4] shows the performance of 
different orders of the power method compared to peeling polynomials of a similar order. 
Clearly the peeling polynomials allow for a more efficient calculation of the eigenvalues in 
the considered region. The improvement of the Arnoldi extraction by the polynomials seems 
to be saturated at a certain order. The extra number of matrix vector multiplications com- 
pensates the focusing and acceleration effect. This saturation happens at higher orders for 
the peeling polynomials than for the power polynomials. In case of the peeling polynomials, 
the saturation depends on the eigenvalue density, since for a larger lattice size it happens 
at a larger order. At a smaller lattice size the eigenvalue density seems to be too low to 
profit from the better focused calculation. Eventually the performance is limited when the 
next region of a high eigenvalue density is reached. In the considered spectra these regions 
form a regular pattern similar to the free theory (cf. Fig. [T]). Therefore, the limiting high 
eigenvalue density can be attributed to the first doublers. A steep rise of computation time 
is visible at this point, especially on larger lattices. 

Besides the time of the computation the required memory can be a limitation of the 
eigenvalue computations. In that respect, the peeling approach exhibits decisive advantages 
with respect to the power method. This is shown in the right part of Fig. [5| where the 
needed memory is represented by the number of vectors used in the computation. With 
respect to this requirement even quite large orders of the peeling polynomials can lead to 
an improvement. 

Using larger orders of the peeling polynomials on smaller lattices one observes that most 

of the time for the Arnoldi computation is spent on a certain set of configurations with 

( f) 

a low eigenvalue density inside ^i^w- avoid this effect one can imply a small limit on 
the maximal number of iterations such that a smaller number of eigenvalues is extracted 
on these rather uninteresting configurations. The resulting improvement is shown in the 
right part of Fig. [5| Nevertheless, one should be careful with the limit on the maximum 
number of iterations. The chance of missing some eigenvalues is increased when the number 
of iterations gets very low. The Arnoldi algorithm requires a balance of the number of 
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FIG. 4. The performance of the different polynomials with different sizes of the eigenspace and 
different maximal number of iterations. The numbers in front of the polynomial name indicate the 
polynomial order. The computations were done on 91 / 11 configurations of one-flavor QCD for 
the 8^ X 16 / 12^ x 16 lattice. All eigenvalues counted for the number of real A have an imaginary 
part smaller than 0.05. The large variation for the order 64 peeling polynomial on the small lattice 
shows a dependence on the maximal number of iterations as in Fig. |5] 



multiplications in the polynomial and of the Arnoldi iterations. 



VI. TECHNICAL DETAILS OF THE POLYNOMIAL DESIGN 

The polynomials can be designed to improve the focusing on the largest or the smallest 
real eigenvalues. To simplify the notation, we assume that in the latter case the transfor- 
mation Dw — i- 2 — Dw is applied, such that again the largest real eigenvalues should be 
computed,^ The region of the wanted eigenvalues is hence ^^^^^ = {Aj| |3 [Aj] | < e, 3? [Aj] > 
^^min}; whcrc £ > Is Small and Xmm is the deepest point in the spectrum considered 
in the computation. For simplicity, the normalization of the polynomial is chosen such 
that p(xmin) = 1 in each step of the peeling transformation. Thus, cTfc is replaced using 
(To = 1 - x^in/ro or cxfc = 1 - l/r^ for A; > 0. 

To understand the effect of one step of the peeling transformation (i. e. a power trans- 
formation), we represent the complex eigenvalues Aj by their radius and phase, 

= ^/(3^^ [A,] /ro - (To)2 + (53 [A,] /ro)^ and 6, = arctan [(3? [A,] - aor^)/'^ [A,]] , 

after the shift and rescaling. The phase is mapped onto ng^j and a fraction of the eigenvalues 
with a nonzero imaginary part are hence "rotated away" from the real axis and out of the 
region of the computed eigenvalues. This effect focuses the calculation on the real eigen- 
modes. However, the eigenvalues with the largest 6 can be "rotated" inside the computed 

^ For D\Y, this is a symmetry transformation and the spectrum remains unchanged, but we applied the 
method also for other operators, like the even-odd preconditioned Wilson-Dirac operator. 
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FIG. 5. The figure on the left hand side illustrates the memory (number of vectors) required for the 
computation of a number of wanted eigenvalues (12^ x 16 lattice). The eigenspace size was set to 
twice the number of considered eigenvalues plus two. On the right, we show the dependence on the 
maximal number of iterations allowed in the Arnoldi algorithm for an order 256 peeling polynomial 
(8^ X 16 lattice). In this figure the number of considered eigenvalues is enlarged, but the maximum 
number of iterations is kept at a small value. 



region. Let ^max be the maximal phase of all Aj with pi > 1. One way to avoid such an 
entering of the "rotated" eigenvalues is to apply the restriction no^'max < 37r/2. 

The focusing effect can be better controlled when it is visualized by a plot of the contour 
3? [p(A)] = 1 in the complex plane. The eigenvalues in the region of all A with 3? [p(A)] > 1, 
i. e. inside the contour, are computed by the algorithm when it reaches the real eigenvalue 
A = Xjnin- There are no of such regions, and the contours surrounding them tend for p — )■ ±oo 
to the lines 6 = (2/+l)7r/ {2no), with / = 0, 1, . . . , no — 1. The larger the number of eigenvalues 
in ^i)^ divided by the number of eigenvalues in these regions, the better is the focusing of 
the polynomial. The restriction to avoid an entering of the "rotated" points in the computed 
region can now be made more precise: the parameters are restricted such that only the 
region of one contour surrounding has overlap with ^^^^^ • This region should be made 
as small as possible for the best focusing. Thus, for a given no, tq must be minimized as 
much as possible without the appearance of a second contour inside 

Increasing no and adjusting tq by this minimization, one observes that the improvement 
of the focusing saturates at larger no- The contour lines become almost parallel equidistant 
lines for large p (see first plot in Fig. [6]). 

The optimization of the focusing can be applied in each step of the peeling strategy: 
one has to choose a power n^ and minimize r^. The polynomial that deviates most from 
the power polynomial of the same order has = 2 for all k. It is shown in the second 
plot of Fig. |6j Compared to the optimal power polynomial of the same order, the region 
{A|3f? [p(A)] > 1} is narrower at the parts of the spectrum that are computed first and slightly 

^ The minimization of tq leads to a maximal 9 and hence to a maximal effect of the "rotation". 
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FIG. 6. The first plot shows the region {A|3? [p(A)] > 1} and the surrounding contours for two 
power polynomials obtained in the optimization procedure (xmin = !)■ The black rectangle shows 
the assumed region ^^^^ and the gray lines are 9 = {21 + l)7r/(2?7-o) with / = 0, 1, . . . , tt-q — 1 for 
no = 12. The second plot shows the result of an iterated optimization of the peeling strategy using 
nfc = 2 for A; = 0, . . . , 5 (overall order of the polynomial 64). The eigenvalues in the regions with 
the darkest green are computed first. 



broader for the inner parts of the spectrum. Keeping the overall order of the polynomial 
(product of the rik) fixed, one can adjust the region inside the contour. A larger no, for 
example, leads to a narrowing of the contour in the inner parts of the spectrum and a 
broadening in the outer parts. For the comparison in Section |V] we have chosen Uk = 2 for 
all k. This seems to be the best choice for the eigenvalues of Dw in the outer part of the 
spectrum. In practice it is profitable to test several different polynomials. 

In practical applications some choices of the polynomial might severely lower the preci- 
sion in the multiplication. It might, therefore, be necessary to adapt the parameters, the 
normalization, and the representation of the polynomial. This problem occurs in particular 
for high orders of the polynomials. 

We have calculated the eigenvalues and eigenvectors of the transformed operators p{Dw) 
using the restarted Arnold! algorithm provided by the ARPACK package [6J . The eigenvalues 
of Dyy were obtained from the eigenvectors. Note that in several calculations we have used 
the even-odd preconditioned Wilson-Dirac operator instead of D^. The eigenvalues of the 
preconditioned operator A- are obtained from the eigenvalues of Dyy using A;^) = 2A,-Af. 
In the region of interest this relation is invertible. 

VII. CONCLUSION 

The polynomials obtained with the peeling strategy lead to an efficient calculation of 
the smallest (or largest) real eigenvalues of the Wilson-Dirac operator with the Arnold! 
algorithm. As we have shown in this work they are better adapted for the eigenvalue 
distribution of this operator than simple power transformations. The efficiency of the peeling 
strategy has two main reasons: it circumvents the saturation of the focusing effect in the 
power transformation and the narrowing of the computed region in the outer parts of the 
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spectrum avoids a calculation of regions with a large eigenvalue density close to the real 
axis. Besides this better focusing effect, it provides also an acceleration of the Arnoldi 
algorithm. We have presented a concrete procedure for the optimization of the parameters 
of the polynomials in Section |VI[ 

We have also tested Faber polynomials |15| for the computation of the lowest eigenvalues. 
They offer an interesting alternative with a similar performance as the peeling polynomials 
in the outer parts of the spectrum. A detailed comparison will be the subject of future work. 

The procedure might be adapted for the eigenvalue distribution of other operators with 
a spectrum in a connected region of the complex plane. 
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